#include <iostream>
#include "spline.h"
#include "Bspline.h"
#include "matplotlibcpp.h"
#include <math.h>

namespace plt = matplotlibcpp;

double func(double x)
{
    return 1 / (1+x*x);
}

double dfunc(double x)
{
    return (-2 * x) / ((1+x*x)*(1+x*x));
}



int main()
{
    vector<double> t1,t2,f1,f2;
    for (int i = 0; i < 10; i++)
    {
        t1.push_back(-5+i);
        f1.push_back(func(t1[i]));
        t2.push_back(-4.5+i);
        f2.push_back(func(t2[i]));
    }
    t1.push_back(5);
    f1.push_back(func(5));
    double dfa = dfunc(-5);
    double dfb = dfunc(5);
    double fa = func(-5),fb = func(5);
    Bspline B1(t1,f1);
    B1.cubic_inter(dfa,dfb);
    Bspline B2(t2,f2);
    B2.quadratic_inter(fa,fb);

//以下为画图部分
    double a = -5, b = 5; 
    int n = 100;    //描点画图的点数
    double figgap = (b - a) / n;
    std::vector<double> u(n+1), y(n+1), y1(n+1), y2(n+1);

    plt::figure();          //创建一张图

    for (int j = 0; j <= n; j++)
    {
        u[j] = a + j * figgap;
        y[j] = func(u[j]);
        y1[j] = B1(u[j]);
        y2[j] = B2(u[j]);
    }  
    plt::plot(u, y,{{"label", "$y = \\frac{1}{1+x^2} $"},{"color","k"},{"ls", "--"}});
    plt::plot(u, y1,{{"label", "cubic Bspline"}, {"ls", "-"}});  
    plt::plot(u, y2,{{"label", "quadratic Bspline"}, {"ls", "-"}});  

        
    plt::xlim(a, b);
    plt::title("Bspline in C"); // set a title
    plt::legend();                // enable the legend

    plt::save("../fig/C-func.png");
    cout << "output : C-func.png" << endl;


    //error
    plt::figure();          //创建一张图

    for (int j = 0; j <= n; j++)
    {
        u[j] = a + j * figgap;
        y1[j] = fabs(B1(u[j]) - func(u[j]));
        y2[j] = fabs(B2(u[j]) - func(u[j]));
    } 
    plt::plot(u, y1,{{"label", "cubic Bspline"}, {"ls", "-"}});  
    plt::plot(u, y2,{{"label", "quadratic Bspline"}, {"ls", "-"}});  
        
    plt::xlim(a, b);
    plt::title("Bspline error in C"); // set a title
    plt::legend();                // enable the legend

    plt::save("../fig/C-error.png");
    cout << "output : C-error.png" << endl;

    //error at exact point
    vector<double> x = {-3.5,-3,-0.5,0,0.5,3,3.5};
    cout << "x = ";
    for (int i = 0; i < x.size(); i++) cout << "\t" << x[i];
    cout << endl;
    cout << "E_c = ";
    for (int i = 0; i < x.size(); i++) cout << "\t" << fabs(B1(x[i]) - func(x[i]));
    cout << endl;
    cout << "E_q = ";
    for (int i = 0; i < x.size(); i++) cout << "\t" << fabs(B2(x[i]) - func(x[i]));
    cout << endl;
    

    std::cout << "finish"  << endl;
}